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1 Introduction 


Over the last few years, considerable progress has been made in the development of tech- 
niques for fabricating high-quality superconducting circuits, and this success, together with 
major advances in the theoretical understanding of quantum detection and mixing at mil- 
limetre and sub millim etre wavelengths [1], has made the development of CAD techniques 
for superconducting nonlinear circuits an important new enterprise. For example, arrays of 
quasioptical mixers are now being manufactured, where the antennas, matching networks, 
filters and superconducting tunnel junctions are all fabricated by depositing niobium and 
a variety of oxides on a single quartz substrate. There are no adjustable tuning elements 
on these integrated circuits, and therefore, one must be able to predict their electrical 
behaviour precisely. This requirement, together with a general interest in the generic be- 
haviour of devices such as direct detectors and harmonic mixers, has lead us to develop 
a range of CAD tools for simulating the large- signed, small-signal, and noise behaviour of 
superconducting tunnel junction circuits. 


2 Large-signal analysis 


To model the behaviour of a quasiparticle mixer, it is first necessary to simulate the large- 
signal steady-state dynamics of the local-oscillator circuit. Once the large-signal operating 
point is known, it is then possible to perturb, either numerically or analytically, the under- 
lying system of equations to gain information about the linear relationships between signal 
and noise variables. 

The main problem is how does one calculate the periodic current that flows through a 
tunnel junction when a periodic voltage is applied? For semiconductor devices this calcula- 
tion is almost always carried out in the time domain, and fast Fourier transforms are used 
to interface the terminal waveforms to the frequency- domain description of the embedding 
circuit. 

Classical resistive mixer diodes are relatively easy to simulate because the induced 
current is an instantaneous function of the terminal voltage. Quantum mixer diodes on 
the other hand are difficult to simulate because the tunnelling current depends on the 
voltage that was across the junction at very long times in the past. In the time domain 
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the current is calculated through an integral which is similar to the convolution integral 
of linear systems theory, and the tunnel junction is characterised by a response function 
which oscillates at the gap frequency with an envelope which decays inversely with time 
at large times. To evaluate the tunnelling current it is necessary to sample the terminal 
voltage at a rate greater than the gap frequency , and to integrate beyond a limit which is 
inversely related to the voltage width of the dc nonlinearity. Time- domain simulations are 
useful for studying the switching behaviour of tunnel junctions, but they are inappropriate 
for studying the steady-state behaviour of RF circuits. 

When a sinusoidal potential is applied to a superconducting tunnel junction, the wave- 
functions associated with the quasiparticle states on one side of the barrier are coherently 
phase modulated. The spectrum of the phase factor is a comb of delta functions whose 
coefficients are the elements of a Bessel-function sequence. The Bessel functions have the 
same argument, determined by the voltage drive level, and consecutive orders ranging from 
some large positive integer to the same negative integer. The trick is to recognize that when 
a periodic potential is applied, the spectrum of the overall phase factor is the convolution 
of the spectra associated with the individual harmonic contributions. Once the spectrum 
of the overall phase factor is known, it is possible to calculate the harmonic phasors of the 
tunnelling current from the dc I-V curve and its Hilbert transform [2]. 

The above procedure describes a way of calculating the periodic current that flows in a 
tunnel junction when a periodic potential is applied. In a real circuit the tunnel junction 
is embedded in a linear network and the problem of determining the various voltages and 
currents is complex. Applying the method of harmonic balance [3] to a generic circuit 
comprising a tunnel junction and a Thevenin voltage source, leads to a system of coupled 
nonlinear algebraic equations. Mathematically, the problem then consists of finding the 
roots of these equations; electrically, the problem is equivalent to searching for a waveform 
that simultaneously satisfies the circuit equations at every harmonic frequency. In nonlinear 
CAD terminology, the scheme is a frequency- domain spectral-balance method, however, 
unlike other versions, the spectral decomposition is based on device physics, rather than on 
expanding the terminal behaviour in a set of basis functions. 

The set of algebraic equations that results from applying the method of harmonic bal- 
ance to a tunnel-junction circuit must, in general, be solved numerically. By repeatedly 
analyzing, in different ways, a large pseudo-random set of tunnel-junction circuits, we have 
investigated the speeds and stabilities of a range of iterative root-finding techniques. A 
comparison of the techniques is shown in Fig. 1, where we have plotted the percentage 
of circuits that converge, and the mean number of iterations taken, as a function of the 
damping factor. The damping factor is a coefficient between 0 and 1 which determines the 
degree to which the result of an iteration influences the next guess. A small value improves 
stability at the expense of reducing the rate of convergence. The solid and dashed lines in 
Fig. 1 correspond to two different quality characteristics. It should be appreciated that the 
plots represent a total of around 20,000 circuit simulations. 

The tunnel junction is a nonlinear admittance in the sense that it is most easy to 
calculate the current in terms of the terminal voltage. However, fixed-point voltage-update 
methods [4] are inappropriate for analyzing tunnel-junction circuits, especially in a common- 
user environment, because they fail to converge when the large- signal harmonic admittances 
of the tunnel junction are much greater than those of the embedding circuit. This problem 
is clearly demonstrated in Fig. 1 where it is seen that the routine will only converge if 
the system is heavily damped. A slightly more sophisticated way of finding the roots is 
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Figure 1: Comparison of various techniques for calculating the large-signal quantum be- 
haviour of superconducting tunnel-junction circuits. 


to use a multi-dimensional variant of the secant method [5]. This method is similar to 
the fixed-point method, in the sense that it is only necessary to calculate the tunnelling 
current once per iteration, however, because coarse derivative information is included one 
might expect the routine to behave more reasonably. Somewhat surprisingly, the routine is 
significantly worse despite the additional information. The problem is caused by the fact 
that, effectively, only the terms on the leading diagonal of the Jacobian matrix are included, 
and coupling between harmonics relies on the current calculations. As long as the current 
at. a given harmonic is most strongly influenced by the voltage at the same harmonic then 
the routine will work well. In a highly nonlinear tunnel-junction circuit, however, there 
is strong coupling between harmonics and the routine is inadequate. Fig. 1 shows the 
behaviour of a harmonic-Newton [6] [7] scheme where the full Jacobian matrix is used. It is 
possible to calculate the Jacobian matrix analytically, however, we prefer to calculate the 
Jacobian matrix using finite differences. Harmonic Newton results in the least number of 
failures: in fact, it finds a solution for 75 % of the circuits studied, and to a large extent 
the st ability of the method is independent of the quality of the junction being investigated. 
The fact that the convergence parameter has little effect on this fraction, together with the 
almost reciprocal dependence of the mean number of iterations, shows that if it is possible 
for the method to find a solution then it will eventually do so. Reducing the damping factor 
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simply reduces the size of the voltage steps taken at each iteration, however, these steps 
are usually in the correct direction. 

We have now performed a very large number of real circuit simulations, and despite 
the fact that 25 % of the randomly generated circuits failed to converge, we have never 
come across a real circuit that has not converged. We have investigated this problem in 
some detail, and we have found that many of the circuits that do not converge are behaving 
in a non-periodic manner. This behaviour usually requires that the embedding circuit 
impedances sire very much larger than the normal-state resistance. 

An alternative approach to finding a root in many dimensions is to recast the problem 
into a multidimensionsil optimization. To do this change, the error function is used to 
construct a scalar quantity that has a global minimum at the required root. One of the 
attractions of optimization is that uninteresting variables, such sis the locsil oscillator drive 
level, csin be eliminated from the smalysis by msiking the variable psirt of the objective 
function. In general, we have found that optimization methods sire slow and should not be 
used unless there is a particulsir reason to do so. Unfortunately, there is insufficient time 
to discuss this more advanced topic in this short paper. 

The results of a typical large-signsil smsdysis sire shown in Fig. 2. The sequence of plots 
shows how the pumped dc I-V curve of a typical Nb-AlOx-Nb tunnel junction evolves as 
the uiCR product is changed. As the capacitance decreases the subgap current increases, 
sind non-classicsd negative differentisd resistance is induced on a number of photon steps. 
Notice that large capacitances sire required before the chsiracteristic relaxes to its constsint 
sinusoidsd-voltage-driven form. Also shown, for comparison, is an smalysis where the har- 
monic feedback is turned off. Curiously, it seems as if internsil harmonic pumping can 
enhance the negative differential resistsmce induced on high-order photon steps — for this 
reason it is possible, in certain circumstsmces, for the small-signal behaviour of a mixer to 
be very sensitive to harmonic impedance levels. 


3 Small-signal analysis 


Once the large-signal behaviour of a mixer has been established, it is possible to calcu- 
late the small-signal and noise performance. The admittance and noise-current correlation 
matrices are determined, in the usual way, through quantum-mechanical generalizations of 
commonly used classical concepts. We then use a selection of linear transforms to reduce 
the admittance parameters to two-port impedance and scattering parameters, and the cur- 
rent. correlation matrix to a noise-temperature matrix from which the standard two-port 
noise parameters can be deduced [8] [9]. It transpires that the whole scheme, both signal 
and noise, can be very elegantly normalized to the gap voltage and gap current of the 
tunnel barrier. The advantages of our generalized approach are that one does not have to 
specify before hand which ports are to be used for the input and output, and one can easily 
calculate the two-port small-signal and noise parameters which can then be loaded into pro- 
prietary microwave circuit simulators for further analysis. For example, we are interested in 
designing mixers that have the first stage of low-noise IF amplification in the mixer block. 
A further advantage of our scheme is that the noise performance is described in terms of 
correlated travelling noise waves, and this approach is an elegant way of considering a mixer 
as an integral part of a quasioptical system along which noise waves propagate. 
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Figure 2: The dc I-V curve of a pumped Nb-AlOx-Nb tunnel junction for different values 
of u>C R product. The first plot does not include internal harmonic pumping. 


4 Mixer simulations 


To date we have studied the large-signal, small-signal, and noise behaviour of mixers by 
adopting the design procedure suggested by Kerr [10]. That is to say, the mixer is operated 
in a double-sideband mode, and the source and load impedances are assumed to be real. The 
ratio of the source and load impedances, which in practice is determined by the geometry 
of the mount, is set at some fixed value. In general, we use a value of unity as a higher 
value tends to degrade the input return loss of the mixer. It is interesting to note, however, 
that it may be possible to choose the ratio so as to minimize the sensitivity of the gain 
to variations in the source resistance. The free parameter, as far as the design process 
is concerned, is the normal-state resistance. Although we assume that the capacitance of 
the tunnel junction is tuned out at the fundamental, we assume that the impedance at 
the harmonics is given by the capacitance of the tunnel junction alone. The effects of 
junction capacitance are considered in companion paper [11], here we simply demonstrate 
the procedure by plotting, in Fig. 3, the transducer gain, noise temperature, input return 
loss, and normalized output impedance of a typical Nb-AlOx-Nb mixer as a function of the 
normalized source resistance; the various curves are for different normalized frequencies ( 
normalized to the gap frequency). It is interesting to note that the overall performance is 
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Figure 3: The transducer gain, noise temperature, input return loss, and normalized output 
impedance of a typical Nb- AlOx-Nb mixer as a function of the normalized source resistance. 


very poor for low values of source resistance, and this is probably the single most important 
reason why mixer performances improve significantly when integrated tuning elements are 
used. 


A useful normalized expression can easily be derived for the source resistance at which 
unity gain, good input match, and minimum noise temperature can be achieved. The 
expression is 


^ = 0 . 5 V --° 92 
R, PH 


78 


V g {mV) 


- 0.92 


( 1 ) 


and it. applies for frequencies between 0.2 and 0.8 of the gap frequency. The exponent is 
slightly different from that given by Kerr and Pan, because we have taken into account the 
fact, that the optimum bias point does not remain in the middle of the first photon step 
below the gap for frequencies greater than about 0.5 of the gap frequency. 

It is well known that if one plots the conversion gain, at a given frequency, as a function 
of the uiC R product, at some point the conversion gain becomes depressed. This can be 
regarded as the frequency at which harmonic effects become significant; or equivalently, the 
frequency at which the five-port, rather than the three-port model should be used. Using 
the above value for the source resistance, we have investigated this behaviour and generated 
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Figure 4: Critical current density against frequency for different values of u>C R product. 
The dotted line shows the optimum current density as a function of frequency. 


the following expression for the optimum u>CR product 


wC R 


opt 


= V PH 75 = 61 


' Vg(mV) 

J(GHz) \ 


^ 0.75 


( 2 ) 


and this in turn generates the following expression for the optimum value of the critical 
current density: J c (Acm~ 2 ) - OAf(GHz) 1 ' 5 . We have assumed / C R„ = 1.8ml' and a 
specific capacitance of 45 f F [im~ 2 . Once again the above requirement is less severe than 
that published by Kerr. Finally, in Fig. 4 we plot the optimum critical current density as 
a function of frequency, and we show lines of constant u jCR product. Above the dotted 
line, harmonic effects are important, whereas below the dotted line harmonic effects can be 
ignored. 
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